TP modèles linéaires - Partie 1
4modIA - 2023-2024
L’objectif de ce TP est d’illustrer les notions abordées dans le chapitre de régression linéaire.
Les librairies R nécessaires pour ce TP :
library(ellipse)
library(leaps)
library(MASS)
library(corrplot)
library(glmnet)
library(coefplot)
library(ggplot2)
library(gridExtra)
library(ggfortify)
library(plotly)
library(reshape2)1 Introduction
La pollution de l’air constitue actuellement une des préoccupations majeures de santé publique. De nombreuses études épidémiologiques ont permis de mettre en évidence l’influence sur la santé de certains composés chimiques comme le dioxyde souffre (SO2), le dioxyde d’azote (NO2), l’ozone (O3) ou des particules en suspension. Des associations de surveillance de la qualité de l’air (Air Breizh en Bretagne depuis 1994) existent sur tout le territoire métropolitain et mesurent la concentration des polluants. Elles enregistrent également les conditions météorologiques comme la température, la nébulosité, le vent, les chutes de pluie en relation avec les services de Météo France… L’une des missions de ces associations est de construire des modèles de prévision de la concentration en ozone du lendemain à partir des données disponibles du jour : observations et prévisions de Météo France. Plus précisément, il s’agit d’anticiper l’occurrence ou non d’un dépassement légal du pic d’ozone (180 \(\mu\)gr/m3) le lendemain afin d’aider les services préfectoraux à prendre les décisions nécessaires de prévention : confinement des personnes à risque, limitation du trafic routier. Plus modestement, l’objectif de cette étude est de mettre en évidence l’influence de certains paramètres sur la concentration d’ozone (en \(\mu\)gr/m3) et différentes variables observées ou leur prévision. Les 112 données étudiées ont été recueillies à Rennes durant l’été 2001.
Les 13 variables observées sont :
- maxO3 : Maximum de concentration d’ozone observé sur la journée en \(\mu\)gr/m3
- T9, T12, T15 : Température observée à 9, 12 et 15h
- Ne9, Ne12, Ne15 : Nébulosité observée à 9, 12 et 15h
- Vx9, Vx12, Vx15 : Composante E-O du vent à 9, 12 et 15h
- maxO3v : Teneur maximum en ozone observée la veille
- vent : orientation du vent à 12h
- pluie : occurrence ou non de précipitations
Les données sont disponibles sur la page moodle du cours. Pour les importer, vous pouvez utiliser la commande suivante :
ozone = read.table("Ozone.txt")Afin de vous familiariser avec les données, faites dans un premier
temps une analyse de statistique descriptive. Vous pouvez utiliser les
fonctions summary(), boxplot(),
pairs(), barplot(), corrplot(),
….
# A COMPLETER - Faire des stat descriptives des données
summary(ozone) maxO3 T9 T12 T15
Min. : 42.00 Min. :11.30 Min. :14.00 Min. :14.90
1st Qu.: 70.75 1st Qu.:16.20 1st Qu.:18.60 1st Qu.:19.27
Median : 81.50 Median :17.80 Median :20.55 Median :22.05
Mean : 90.30 Mean :18.36 Mean :21.53 Mean :22.63
3rd Qu.:106.00 3rd Qu.:19.93 3rd Qu.:23.55 3rd Qu.:25.40
Ne9 Ne12 Ne15 Vx9
Min. :0.000 Min. :0.000 Min. :0.00 Min. :-7.8785
1st Qu.:3.000 1st Qu.:4.000 1st Qu.:3.00 1st Qu.:-3.2765
Median :6.000 Median :5.000 Median :5.00 Median :-0.8660
Mean :4.929 Mean :5.018 Mean :4.83 Mean :-1.2143
3rd Qu.:7.000 3rd Qu.:7.000 3rd Qu.:7.00 3rd Qu.: 0.6946
Vx12 Vx15 maxO3v vent
Min. :-7.878 Min. :-9.000 Min. : 42.00 Length:112
1st Qu.:-3.565 1st Qu.:-3.939 1st Qu.: 71.00 Class :character
Median :-1.879 Median :-1.550 Median : 82.50 Mode :character
Mean :-1.611 Mean :-1.691 Mean : 90.57
3rd Qu.: 0.000 3rd Qu.: 0.000 3rd Qu.:106.00
pluie
Length:112
Class :character
Mode :character
[ getOption("max.print") est atteint -- ligne 1 omise ]2 Régression linéaire simple
Dans cette section, nous souhaitons expliquer la concentration d’ozone maximale de la journée (maxO3) en fonction de la température à 12h (T12).
- Représentez le nuage de points \((x_i,y_i)\) à l’aide de la fonction
plot()(ougeom_point()de ggplot2).
# A completer
# Scatter plot
plot(ozone$T12, ozone$maxO3,
xlab = "Temperature at 12:00 PM (T12)",
ylab = "Maximum Ozone Concentration (maxO3)",
main = "Scatter Plot of maxO3 vs. T12")
Question : Pensez-vous que l’ajustement d’un modèle de régression linéaire \(y_i = \theta_0+\theta_1 x_i +\varepsilon_i\) est justifié graphiquement?
Réponse : On voit une espèce de tendance linéaire : un modèle de régression linéaire est justifiée.
- Effectuez la régression linéaire à l’aide de la fonction
lm()et exploitez les résultats.
# Régression linéaire simple
reg.simple = lm(maxO3 ~ T12, data = ozone)
# Afficher un résumé des résultats de la régression
summary(reg.simple)
Call:
lm(formula = maxO3 ~ T12, data = ozone)
Residuals:
Min 1Q Median 3Q Max
-38.079 -12.735 0.257 11.003 44.671
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -27.4196 9.0335 -3.035 0.003 **
T12 5.4687 0.4125 13.258 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 17.57 on 110 degrees of freedom
Multiple R-squared: 0.6151, Adjusted R-squared: 0.6116
F-statistic: 175.8 on 1 and 110 DF, p-value: < 2.2e-16Interprétation :
Que contiennent les sorties suivantes :
reg.simple$coefficients
reg.simple$residualsRéponse : coefficients -> vecteur des valeurs observées de \(\^{\theta}^{obs}\) residuals -> vecteur des résidus
- A l’aide des commandes suivantes, tracez l’estimation de la droite de régression sur le nuage de points ainsi qu’un intervalle de confiance à \(95\%\) de celle-ci :
ggplot(ozone, aes(T12, maxO3))+
geom_point() +
geom_smooth(method=lm, se=TRUE)+
xlab("T12")+ ylab("maxO3")
Ce graphique permet visuellement de vérifier l’ajustement des données au modèle de régression linéaire. Que remarquez-vous?
- A l’aide de la commande suivante, étudiez les résidus
autoplot(reg.simple,which=c(1,2,4),label.size=2) 
On prendra soin de comprendre les différents graphiques obtenus.
- On s’intéresse maintenant à la qualité de prédiction du modèle. On va donc tracer un intervalle de confiance des prédictions à \(95\%\) avec les commandes suivantes. Commentez.
temp_var = predict(reg.simple,
interval="prediction")
new_df = cbind(ozone, temp_var)
ggplot(new_df, aes(T12, maxO3))+
geom_point() +
geom_line(aes(y=lwr), color = "red",
linetype = "dashed")+
geom_line(aes(y=upr), color = "red",
linetype = "dashed")+
geom_smooth(method=lm, se=TRUE)+
xlab("T12")+
ylab("maxO3")- On va maintenant s’intéresser à la construction d’intervalles de
confiance pour \(\theta_0\) et \(\theta_1\) à \(95\%\).
A l’aide de la fonctionconfint(), construisez un intervalle de confiance pour chaque paramètre séparément.
# A COMPLETER
# Intervalles de confiance pour theta_0 (intercept) et theta_1 (pente)
IC = confint(reg.simple)
ICPour tenir compte de la dépendance entre \(\theta_0\) et \(\theta_1\), on peut aussi construire une région de confiance pour le vecteur \(\theta=(\theta_0,\theta_1)'\) grâce aux commandes suivantes. Comparez les résultats.
df1 = as.data.frame(
rbind(coefficients(reg.simple),
ellipse(reg.simple,level=0.95)))
colnames(df1) = c("intp", "slope")
ggplot(data=df1[-1,],aes(x=intp, y=slope))+
geom_path()+
annotate("rect",xmin=IC[1,1],xmax=IC[1,2],
ymin=IC[2,1],ymax=IC[2,2],fill="red",alpha=0.1)+
geom_point(data=df1[1,], aes(x=intp, y=slope),
pch=3)3 Régression linéaire multiple
Dans cette partie, nous souhaitons analyser la relation entre le maximum journalier de la concentration d’ozone (maxO3) et
- la température à différentes heures de la journée (T9, T12, T15),
- la nébulosité à différentes heures de la journée (Ne9, Ne12, Ne15),
- la projection du vent sur l’axe Est-Ouest à différentes heures de la journée (Vx9, Vx12,Vx15),
- la concentration maximale d’ozone de la veille (maxO3v)
On va donc utiliser le sous-jeu de données suivant
ozone1 = ozone[,1:11]et mettre en place un modèle de régression linéaire multiple.
- Faites une analyse descriptive de ce jeu de données.
# Afficher un résumé statistique des variables
summary(ozone1)
# Afficher des boîtes à moustaches pour visualiser la distribution des variables numériques
boxplot(ozone1)
# Afficher une matrice de nuages de points pour explorer les relations entre les variables
pairs(ozone1)- Rappelez l’écriture mathématique du modèle de régression linéaire multiple.
\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + ... + \beta_p X_p + \varepsilon \]
- Ajustez un modèle de régression linéaire multiple à l’aide de la
commande
lm()et commentez les résultats (on appelle reg.mul la sortie de R dans la suite).
reg.mul = lm(maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v, data = ozone1)
resume.mul = summary(reg.mul)
summary(reg.mul)
Call:
lm(formula = maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 +
Vx12 + Vx15 + maxO3v, data = ozone1)
Residuals:
Min 1Q Median 3Q Max
-53.566 -8.727 -0.403 7.599 39.458
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 12.24442 13.47190 0.909 0.3656
T9 -0.01901 1.12515 -0.017 0.9866
T12 2.22115 1.43294 1.550 0.1243
T15 0.55853 1.14464 0.488 0.6266
Ne9 -2.18909 0.93824 -2.333 0.0216 *
Ne12 -0.42102 1.36766 -0.308 0.7588
Ne15 0.18373 1.00279 0.183 0.8550
Vx9 0.94791 0.91228 1.039 0.3013
Vx12 0.03120 1.05523 0.030 0.9765
Vx15 0.41859 0.91568 0.457 0.6486
maxO3v 0.35198 0.06289 5.597 1.88e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 14.36 on 101 degrees of freedom
Multiple R-squared: 0.7638, Adjusted R-squared: 0.7405
F-statistic: 32.67 on 10 and 101 DF, p-value: < 2.2e-16- Etudiez les résidus
autoplot(reg.mul,which=c(1,2,4),label.size=2) 
4 Sélection de variables
Dans le summary(reg.mul), un test est fait sur chaque
coefficient. Il revient à tester que la variable n’apporte pas
d’information supplémentaire sachant que toutes les autres variables
sont dans le modèle. Il est donc préférable d’utiliser des procédures de
choix de modèles pour sélectionner les variables significatives. On va
ici comparer la sélection de variable obtenue par différents critères:
BIC, AIC, \(R^2\) ajusté, Cp de
Mallows. Pour cela, vous pouvez utiliser la fonction
regsubsets() de la librairie leaps et la fonction
stepAIC() de la librairie MASS. Commentez les
résultats obtenus avec les différents critères, vous pourrez vous aider
des commandes suivantes :
# Utiliser regsubsets pour sélectionner les variables avec différents critères
choix = regsubsets(maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v, data = ozone1)
plot(choix, scale = 'Cp')
plot(choix, scale = 'adjr2')
plot(choix, scale = 'bic')
# Utiliser stepAIC pour sélectionner le modèle basé sur l'AIC
stepAIC(reg.mul)Start: AIC=607.26
maxO3 ~ T9 + T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 +
maxO3v
Df Sum of Sq RSS AIC
- T9 1 0.1 20827 605.26
- Vx12 1 0.2 20827 605.26
- Ne15 1 6.9 20834 605.30
- Ne12 1 19.5 20847 605.36
- Vx15 1 43.1 20870 605.49
- T15 1 49.1 20876 605.52
- Vx9 1 222.6 21050 606.45
<none> 20827 607.26
- T12 1 495.5 21323 607.89
- Ne9 1 1122.6 21950 611.14
- maxO3v 1 6459.6 27287 635.51
Step: AIC=605.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx12 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Vx12 1 0.1 20827 603.26
- Ne15 1 6.9 20834 603.30
- Ne12 1 22.1 20849 603.38
- Vx15 1 43.5 20871 603.49
- T15 1 49.4 20877 603.52
- Vx9 1 264.0 21091 604.67
<none> 20827 605.26
- T12 1 641.4 21469 606.66
- Ne9 1 1183.6 22011 609.45
- maxO3v 1 7112.2 27940 636.16
Step: AIC=603.26
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Ne15 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Ne15 1 6.9 20834 601.30
- Ne12 1 23.1 20850 601.38
- T15 1 50.0 20877 601.53
- Vx15 1 79.6 20907 601.69
- Vx9 1 320.1 21148 602.97
<none> 20827 603.26
- T12 1 647.3 21475 604.69
- Ne9 1 1185.1 22013 607.46
- maxO3v 1 7112.6 27940 634.16
Step: AIC=601.3
maxO3 ~ T12 + T15 + Ne9 + Ne12 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Ne12 1 16.2 20851 599.38
- T15 1 45.2 20880 599.54
- Vx15 1 75.1 20910 599.70
- Vx9 1 325.2 21160 601.03
<none> 20834 601.30
- T12 1 971.8 21806 604.40
- Ne9 1 1215.8 22050 605.65
- maxO3v 1 7106.3 27941 632.17
Step: AIC=599.38
maxO3 ~ T12 + T15 + Ne9 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- T15 1 45.6 20896 597.63
- Vx15 1 77.6 20928 597.80
- Vx9 1 337.6 21188 599.18
<none> 20851 599.38
- T12 1 1034.0 21885 602.80
- Ne9 1 2155.5 23006 608.40
- maxO3v 1 7130.2 27981 630.33
Step: AIC=597.63
maxO3 ~ T12 + Ne9 + Vx9 + Vx15 + maxO3v
Df Sum of Sq RSS AIC
- Vx15 1 74.1 20970 596.02
- Vx9 1 366.5 21263 597.58
<none> 20896 597.63
- Ne9 1 2241.1 23137 607.04
- T12 1 6712.6 27609 626.83
- maxO3v 1 7381.5 28278 629.51
Step: AIC=596.02
maxO3 ~ T12 + Ne9 + Vx9 + maxO3v
Df Sum of Sq RSS AIC
<none> 20970 596.02
- Vx9 1 903.4 21874 598.75
- Ne9 1 2714.8 23685 607.66
- T12 1 6650.4 27621 624.88
- maxO3v 1 7363.5 28334 627.73
Call:
lm(formula = maxO3 ~ T12 + Ne9 + Vx9 + maxO3v, data = ozone1)
Coefficients:
(Intercept) T12 Ne9 Vx9 maxO3v
12.6313 2.7641 -2.5154 1.2929 0.3548 Avec le critère BIC, nous retenons 4 variables : T12, Ne9, Vx9 et
maxO3v.
A l’aide des commandes suivantes, testez le sous-modèle avec les 4
variables retenues par BIC contre le modèle complet. Commentez.
reg.fin=lm(maxO3~T12+Ne9+Vx9+maxO3v,data=ozone1)
summary(reg.fin)
anova(reg.fin,reg.mul)5 Régressions régularisées
On commence par centrer et réduire les données avant de mettre en place et comparer des méthodes de régression régularisées.
tildeY=scale(ozone1[,1],center=T,scale=T)
tildeX=scale(ozone1[,-1],center=T,scale=T)5.1 Régression Ridge
- A l’aide de la fonction
glmnet(), ajustez une régression ridge en faisant varier \(\lambda\) sur une grille. On stockera le résultat dans la variablefitridge. Explorez le contenu defitridge.
lambda_seq<-10^(seq(-4,4,0.01))
# Ajuster la régression ridge avec glmnet
fitridge <- glmnet(x = tildeX, y = tildeY, alpha = 0, lambda = lambda_seq, family = c("gaussian"), intercept = F)
summary(fitridge) Length Class Mode
a0 801 -none- numeric
beta 8010 dgCMatrix S4
df 801 -none- numeric
dim 2 -none- numeric
lambda 801 -none- numeric
dev.ratio 801 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 7 -none- call
nobs 1 -none- numeric- Tracez les chemins de régularisation de chaque variable et commentez.
df=data.frame(lambda = rep(fitridge$lambda,ncol(tildeX)), theta=as.vector(t(fitridge$beta)),variable=rep(colnames(tildeX),each=length(fitridge$lambda)))
g1 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
geom_line()+
theme(legend.position="bottom")+
scale_x_log10()
ggplotly(g1)- A l’aide de la fonction
cv.glmnet()mettez en place une validation croisée pour sélectionner le “meilleur” \(\lambda\) par MSE.
# Réaliser une validation croisée pour sélectionner le meilleur lambda
ridge_cv <- cv.glmnet(x = tildeX, y = tildeY, alpha = 0, lambda = lambda_seq, nfolds = 10, type.measure = c("mse"), intercept = F)
# Trouver la meilleure valeur de lambda en minimisant le MSE
best_lambda <- ridge_cv$lambda.mindf2=data.frame(lambda=ridge_cv$lambda,MSE=ridge_cv$cvm,cvup=ridge_cv$cvup,cvlo=ridge_cv$cvlo)
gmse<-ggplot(df2)+
geom_line(aes(x=lambda,y=MSE))+
geom_vline(xintercept = ridge_cv$lambda.min,col="red",linetype="dotted")+
geom_line(aes(x=lambda,y=cvup),col="blue",linetype="dotted")+
geom_line(aes(x=lambda,y=cvlo),col="blue",linetype="dotted")+
#xlim(c(0,ridge_cv$lambda.min+0.5))+
scale_x_log10()
ggplotly(gmse)La valeur de \(\lambda\) sélectionnée vaut
g1=g1 +
geom_vline(xintercept = best_lambda,linetype="dotted", color = "red")+
scale_x_log10()5.2 Régression Lasso
- A l’aide de la fonction
glmnet(), ajustez une régression Lasso en faisant varier \(\lambda\) sur une grille. On stockera le résultat dans la variablefitlasso. Explorez le contenu defitlasso.
# Ajuster la régression Lasso avec glmnet
fitlasso <- glmnet(x = tildeX, y = tildeY, alpha = , lambda = lambda_seq, family = c("gaussian"), intercept = F)
# Afficher un résumé des résultats
summary(fitlasso) Length Class Mode
a0 801 -none- numeric
beta 8010 dgCMatrix S4
df 801 -none- numeric
dim 2 -none- numeric
lambda 801 -none- numeric
dev.ratio 801 -none- numeric
nulldev 1 -none- numeric
npasses 1 -none- numeric
jerr 1 -none- numeric
offset 1 -none- logical
call 6 -none- call
nobs 1 -none- numeric- Tracez le chemin de régularisation de chacune des variables et commentez
df=data.frame(lambda = rep(fitlasso$lambda,ncol(tildeX)), theta=as.vector(t(fitlasso$beta)),variable=rep(colnames(tildeX),each=length(fitlasso$lambda)))
g3 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
geom_line()+
theme(legend.position="bottom")+
scale_x_log10()
ggplotly(g3)- A l’aide de la fonction
cv.glmnet()mettez en place une validation croisée pour sélectionner le “meilleur” \(\lambda\) par MSE. En pratique, il est préconisé d’utilisélambda.1se(la plus grande valeur de \(\lambda\) telle que l’erreur standard se situe à moins de 1 de celle du minimum).
lasso_cv <- cv.glmnet(...) # A COMPLETER
best_lambda <-lasso_cv$lambda.min
lambda1se <- lasso_cv$lambda.1seLa valeur de \(\lambda\) sélectionnée vaut ….
g3=g3 +
geom_vline(xintercept = best_lambda,linetype="dotted", color = "red")+
geom_vline(xintercept = lambda1se,linetype="dotted", color = "blue")+
scale_x_log10()
g3- Quelle sélection de variables obtient-on alors ?
Vous pouvez utiliser la fonction extract.coef() de la
librairie coefplot
# A COMPLETER5.3 Régression Elastic Net
- Reprenez les questions précédentes pour ajuster une régression Elastic Net
fitEN <- glmnet(...)
df=data.frame(lambda = rep(fitEN$lambda,ncol(tildeX)), theta=as.vector(t(fitEN$beta)),variable=rep(c(colnames(tildeX)),each=length(fitEN$lambda)))
g4 = ggplot(df,aes(x=lambda,y=theta,col=variable))+
geom_line()+
theme(legend.position="bottom")+
scale_x_log10()EN_cv <- cv.glmnet(...) # A COMPLETER
best_lambda <-EN_cv$lambda.min
g4=g4 + geom_vline(xintercept = best_lambda,linetype="dotted",
color = "red")
ggplotly(g4)- Comparez les coefficients obtenus avec la régression linéaire, la régression ridge, la régression Lasso et la régression ElasticNet